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Abstract 

We propose a novel method to describe realistically ionization processes with absorbing bound- 
ary conditions in basis expansion within the formalism of the so-called Non-Adiabatic Quantum 
Molecular Dynamics. This theory couples self-consistently a classical description of the nuclei with 



' a quantum mechanical treatment of the electrons in atomic many-body systems. In this paper we 

o : 

■ extend the formalism by introducing absorbing boundary conditions via an imaginary potential. 

o . 

It is shown how this potential can be constructed in time-dependent density functional theory in 



basis expansion. The approach is first tested on the hydrogen atom and the pre-aligned hydro- 



gen molecular ion H^" in intense laser fields where reference calculations are available. It is then 
Q" 1 applied to study the ionization of non-aligned and H2. Striking differences in the orientation 

dependence between both molecules are found. Surprisingly, enhanced ionization is predicted for 
' perpendicularly aligned molecules. 
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I. INTRODUCTION 



The experimental and theoretical investigation of the interaction of atoms, molecules 
and clusters with intense laser fields represents one of the most challenging problems of 



current research. In atoms, high harmonic generation (HHG 



ionization 



3,0,3 



or stabilization against ionization 



1 3, 3, 



above threshold 



M have been observed. In 



molecules, due to the additional nuclear degrees of freedom (DOF), further mechanisms 




ike molecular stabilization against dissociation 



bond softening 



above threshold dissociation 



22j, orientation dependent HHG 



species dependent orientation dependence of single ionization in dimers |32j, |3 




□ 



14, 



, c harge r esonance enhanced ionization 



molecular 
or an 



unexpected 



atom 13 



s upp ression of ionization in dimers in comparison to the so-called companion 
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38. 



39 



40 



41 



42 



to name but a few effects. 



The theoretical understanding of these mechanisms requires, in principle, the solution 



ectrons and all nuclear DOF. 
, atomic helium 4|| |^ | , 



4JJ45L; 



of the time-dependent Schrodinger equation (TDSE) for al 
However, only for the smallest systems, atomic hydrogen 
laser aligned H^~ j3| and laser aligned H 2 with fixed nuclei ; 
the TDSE exist. In larger systems (or for hydrogen molecules without constraints) approx- 
imations are necessary due to the exponential scaling of the computational effort with the 
number of DOF. 

One possibilty consists in the combination of time-dependent density functional the 



ory |5l| (TDDFT) and a classical description of the nuclei 
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hi 
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56 



57 



58. 



59]. 



Most of the approaches use a representation of the 
solve the time-dependent KS equations 



Kohn-S 



56 



57 J5l 

i. 2 



58, 



lam (KS) functions on a grid to 



approaches has been developed by Reinhard et al. 
and Castro et al. 3 have also developed such methods 



56 



5£ |. The first of these grid based 



. More recently, Dundas 



In contrast, a basis expansion of the KS-orbitals with local basis functions is used in 



the so-ca 
group 



52 



led non-adiabatic quantum molecular dynamics (NA-QMD), developed in our 



531 ] . It has been successfully applied so far to very different non-adiabatic pro- 
cesses, like atom-cluster collisions j^J, ion-fullerene collisions [fnj], laser induced excitation 
and fragmentation of molecules (3 or fragmentation and isomerization of organic molecules 
in laser fields J3- However, a realistic description of ionization with the NA-QMD theory 
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is still an open problem. 

In this work, we present a general method in basis expansion and extend the NA-QMD 
formalism to describe ionization in many-electron systems. To this end, two problems have 
to be addressed. First, an appropriate basis set suitable for the description of highly excited 
and ionized states has to be found. We have focused on this problem in a recent paper . 
The second problem concerns the introduction of absorbing boundary conditions. This is 
still a challenging problem for many-electron calculations performed on grids (see e.g. 
and in particular, completely open for methods using basis expansion. 

In |5J] we have used an ad-hoc manipulation of the electronic expansion coefficients in 
one electron calculations. In this paper, we introduce general absorbing boundary conditions 
in basis expansion via an imaginary potential. It will be shown, how such an potential can 
be build using the many-body Schrodinger equation as well as TDDFT. 

The outline of this paper is as follows. First, we introduce absorbing boundary conditions 
in basis expansion for the general case of the Schrodinger equation in section III Al In sec- 
tion lll B{ the extension of the NA-QMD formalism including absorbing boundary conditions 
is presented. Details of the used absorbing potential are outlined in section Hi CI The method 
is tested on the hydrogen atom (section IIII Aj) as well as aligned H J (section IIIIB|) where 
reference calculations are available. In section IIII CI the calculated ionization probabilities 
and rates of non-aligned and H 2 are presented. A completely different orientation de- 
pendence between both molecules is found. In addition, the calculations predict surprisingly 
enhanced ionization for perpendicularly aligned molecules. 

II. ABSORBING BOUNDARY CONDITIONS 
A. The Schrodinger equation and basic idea 

We start with the introduction of absorbing boundary conditions for the general case of 
the Schrodinger equation (atomic units are used throughout the paper) 




(i) 



where the Hamiltonian 



Hit) 



f + V(t) 



(2) 
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consists of the operator for the kinetic energy t and the potential V which includes the two- 
body interaction. Thus, represents the, in general, many-particle wave function. The 
absorbing boundary conditions are incorporated via an imaginary potential. The Hamilto- 
nian is thus modified, 

H*»{t) = H{t)-iV^{t) (3) 

where Vabs is the absorbing potential. With a hermitian V^bs such an operator if abs is not 
hermitian and therefore the norm is not conserved, i.e. 

A<v^) = (*|i£: bs -i£ abs |*> 



dt 



-2(*|K bs |^). (4) 



A semi-positive definite V^ bs leads to the desired effect of norm reduction, i.e. the derivative 
in (j3J) is negative or zero. 

The balance of the total energy 

E(t) = <*(*) \H(t)Mt)) (5) 



is changed if the absorber potential is used in the propagation of i.e. 

dt v 'dt 



^-E = (*\±vm + A abs (6) 



with 

A abs = -(V\V abs H + HV ahs \V) . (7) 

The additional term is due to the absorber potential and changes the energy balance. Its 
actual effect depends on the definition of the potential 14bs- 
Introducing the time-dependent eigenfunctions \\ a ) to H 

H(t)\ X a(t)) = E a (t)\ Xa (t)) (8) 

we define the absorbing potential as 

oo 

V ahs = ^ fa\Xa) (Xa\ ■ (9) 
a=l 

The states \Xa)i sometimes called "field-following" adiabatic states form an or- 

thonormal set. With our definition of K, bs these states \x a ) are also eigenstates to ff abs 



H ahs (t)\Xa(t)) = (E a (t)-if a )\ Xa (t)) 
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(10) 



but lead to imaginary eigenenergies, i.e. finite lifetimes. The factors f a determine the 
strength of the absorber at a certain energy and are discussed in section III CI The wave 
function \^(t)) is now expanded also in these eigenfunctions 

oo 

l*) = 5> a lxa>. ( n ) 

a=l 

Inserted into the time-derivative of the norm (J3J) this yields 

, oo 

-(^) = -2j2\a a \ 2 fa (12) 

a=l 

and the additional term A abs ((7j) of the energy balance becomes 

oo 

Aabs = -2^\a a \ 2 f a E a . (13) 
0=1 

From equation (|12|). one can see that the absorber potential decreases the norm of arbitrary 
wave functions 1^) only, if all f a > 0. Furthermore, it has to be guaranteed that electronic 
density in bound states is not affected by the absorbing potential. In calculations on spatial 
grids this is approximately satisfied by applying the absorbing boundary conditions far away 
from the nuclei (see e.g. 4jj |6y, |67[). In our case of basis expansion (|11|). this condition can 
naturally be fulfilled if the time- dependent eigenvalues f a of \4bs © are chosen to be 

f/. = if£ a <0 

fa=< , (14) 

(f a >0 if E a >0 

i.e. the absorbing potential acts only on states in the continuum. Thus, A a b s is always zero 
or negative if the potential is defined as in (JHJ) and if the f a meet the criterion (jUj) . 



B. Non-Adiabatic Quantum Molecular Dynamics 

So far we have shown how to introduce absorbing boundary conditions if the Schrodinger 
equation in basis expansion is used. We show in the following how to introduce an absorbing 
potential within the NA-QMD formalism. 

The coupled equations of motion (EOM) for the nuclear and electronic system have been 

n n 

given elsewhere in TDDFT and time-dependent Hartree-Fock (TDHF) |£2j and will not 
be repeated here. Instead we will present the changes in the single-particle EOM arising 
from an additional imaginary potential in the effective single-particle potential. 



The single-particle wave functions are expanded in a local basis 

V°(r, t) = J]<(t)0 a (f- R Aa (t)) with j = 1 . . .N:, a =T, | 



(15) 



a=l 



and only the expansion coefficients a J ^(t) are explicitly time-dependent. The symbol A a 
denotes the atom to which the atomic basis function Q is attached. The Ny, basis functions 
are either located at the nuclei, which in general move, or are located at fixed positions in 
space 6^ . N£ is the number of electrons with the spin a. 

The variation of the total action with respect to electronic expansion coefficients and 
nuclear coordinates leads to self-consistently coupled EOM H; Q] for the electronic ex- 
pansion coefficients a? a a (t) and for the nuclear coordinates Ra- Here we add an imaginary 
potential to the effective single-particle Hamiltonian. The EOM for the electronic expansion 
coefficients are then modified (cf. with 



,J rJ 



dt 



In dini) 



07 



j = l,...,N°,a=1,l ■ 



is the overlap matrix between basis states, 



(16) 



(17) 



Ha/3 - ' 



H, 



cff 



(18) 



is the effective Hamilton matrix with H° s the effective one-particle Hamiltonian 

/ d 

B af3 = ( (p a —(p 

is the non-adiabatic coupling matrix and 



ian 



(19) 



V. 



abs,a/3 



^abs 



(20) 



is the matrix element of the additional absorber potential introduced here and still to be 
defined (see below). 

At this point we note explicitely that the EOM (jlfjjl are exactly the same for TDDFT 
and TDHF Q- The difference between both approaches consists in the calculation of the 
matrix elements H°g which in case of TDHF contains the non-local exchange term. Thus, 



6 



the whole following discussions and derivations belong simultaneously to both approaches, 
TDDFT 3 and TDHF 6^- 

With the additional damping term (J2U)) the time-dependence of the norm (i.e. the total 
number of electrons in this case) becomes 



jo* 3 a 
/3 u a u /3 



(21) 



The total energy reads 



A=l 



E(t) = u(R,t) + j2^R 2 A + E EE<* T ^ a ; 

o-=T, I J =1 a/? 



'/9 



y d 3 rp(f, £) n 



+ /d J rp(f,t)(y(f, J R,t) + - / 



p(f',£) 



£dV) +E xc [p](t) 



(22) 



— # — * 

with the kinetic and potential energy U(R,t) of the nuclei, the external potential V(r, R, t) 
that contains the electron-nuclear interaction and e.g. a laser field Vl(£), the electronic 
density 

*) = E E ¥<T *^ t)^ j<T (f, t) , (23) 



the matrix element of the kinetic energy of the electrons 



T, 



ml 



A 
~2 



(24) 



and the exchange and correlation energy E xc . With the EOM (|16|) and the classical EOM 



for the nuclei 
energy balance 



53j, which are not changed by the imaginary potential, one obtains for the 



-E 



p(f, t) 



dV L (r,t) 
dt 



d 3 r-^Z A 

A=X 



dVL(R A ,t) 
dt 



A 



abs 



(25) 



al a *ai a . (26) 



with Vl the external, time-dependent potential (e.g. a laser) and the additional term 

A abs = - E E E (^bs, a 7 ( 5_1 ) 7 5 H h + ( 5_1 ) 7 5 Kbs, 

The first two terms arise naturally from the interaction of the electronic density p(r, t) and 
the nuclei (with charge Z &) with an external field. They are of course identical with the 
energy balance given in 53j . The last term A a b s in (J25j) is evidently induced by the imaginary 
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potential. We note explicitely that the non-adiabatic coupling matrix B a p (fT9~j). which at 
first glance seems to act equivalently to the absorbing part V^ s a/3 in (fT5j). does not affect 
the energy balance because these terms are canceled out in the calculation of 4? due to the 

OJ I I at 

classical EOM as it should be and as it has been shown in [53j |. 

The general results (|21jl and (|2fij) are equivalent to (j3J) and (jZJ). They are valid for any 
absorbing potential V^ s and any basis set The same holds for the EOM (|16j). 

However, physically, any choice of V^ s must guarantee that the absorption is only applied 
to that part of the density which belongs to the continuum. This is equivalent to the 
requirement that norm and energy are decreased for any arbitrary density (cf. also with 
section III AJ) . i.e. that 

A abs < and ^N<0 V<. (27) 

In order to realize that we make use of the general idea presented in section III Al 

To this end, the single-particle wave functions j^ ') are expanded in the, now, effective 
single-particle "field- following" adiabatic states, i.e. (cf. equation (|11|)) 

N b 

\W){t) = Y J <{t)\ Xa ){t)) with j = l...N:,a=ti ■ (28) 

a=l 

with |Xo(t)) defined as (cf. equation (JBJ) 

HlMxl){t) = el{t)\xl){t). (29) 
The absorber potential is formally constructed as before (cf. equation (JHJ)) 

a=l 

In principle, the expansion coefficients a 3 a a in equation (j28j) can be obtained from the 
EOM ()16|1 written for the basis (|28|). In this case, the basis functions themselves would 
depend on the effective Hamiltonian H° s . Alternatively, and this is done in the present 
work, one may also determine the coefficients a?° by solving (|29J) as a generalized eigen- 
value problem and making use of transformations between the basis sets (J28|) and (fT5|) (see 
appendix EJ • 

In the basis (128)1 and with the absorber ((HP)) the derivative of the norm (j2"Tj) becomes 



^ = - 2 EEEKbs, aa ki 2 (si) 



EE 

8 



and the additional term f|26[) of the energy balance (|25|) is 



A abs = -2 E E ^bs, I «r 1 2 • (32) 

o-=T,J. J=l »=o 

Thus, both quantities are always zero or negative if 

jo for e a a < 

K,bs, aa = fa = \ • (33) 

I > for e£ > 

It now becomes apparent that our choice of the absorber potential does guarantee the phys- 
ical requirement, namely, that density is removed only from states which contribute to the 
continuum. The eigenvalues f a have still to be determined which will be done in the next 
section. 



C. The parameters of the absorber 



The f a are directly connected to the lifetime r a of the states |% a ) if the Hamiltonian H 
is time-independent, i.e. 



fa = + 



1 

2^ 



(34) 



We use here the quantities r a for an appropriate parameterization of f a . It is natural to 
assume that the lifetimes increase smoothly with energy. Thus, we parameterize the lifetimes 
according to 

oo e a < a.u. 

,7 nin \ < e a < E TC f 

T~m'm ^ -E'ref 



> 



(35) 



where the e a are the "field-following" time-dependent energies ()29|) and r min and E vei are two 
parameters still to be determined. 

The general conditions that the parameters r m ; n and E Ie f (or the absorbing potential 
at all) have to satisfy are similar to those in calculations on spatial grids. In this case, 
the absorber has be strong enough to prevent unphysical reflections of the density at the 
boundaries of the grid. On the other side, it must be weak enough to prevent reflections at 
the absorbing potential itself In our case, this is equivalent to the requirements, that 
the absorption is strong enough to prevent reflections at the boundary of the Hilbert space 
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(denned by the finite number of basis states). On the other side, it should be weak enough 
to avoid any suppression of excitation. 

To illustrate this somewhat exceptional, quantum mechanical property one can consider 
a two level system 

i^ai(*) = E 1 {t)a 1 {t) + H 12 (t)a 2 {t) , (36) 

i— a 2 {t) = E 2 {t)a 2 {t) + H 12 (t) ai {t) (37) 

where a\i 2 are the expansion coefficients of the two states, Eu 2 the energies and H\ 2 is the 
coupling matrix element. It is assumed without loss of generality that the basis states 
and |<&2) are orthogonal. The population of the first state | czi | 2 is constant if the absorption 
in the second state is so strong that a 2 = for all times, i.e. 

^ t \ ai (t)\ 2 = iH 12 (t)(a 2 (tya 1 (t)-a 1 (tya 2 (t)) = Q if a 2 (t) = 0W. (38) 

Obviously, in this case, the absorber completely prevents any excitation. 

In the following, we use = 0.3 a.u. because the underlying Hilbert space (i.e. the finite 
basis set used) yields a dense spectrum of states up to this energy. This fulfilles the second 
requirement, discussed above. The minimal decay time, which determines the strength of 
the absorber, is fixed to be r m i n = 5 a.u. In our test calculations we found, however, that 
a relatively large range of parameters leads to very similar results. In addition, the above 
given and fixed parameters are applicable in a very large range of laser parameters as will 
be shown in the following by comparing the results with that of reference calculations. 



III. RESULTS 



A. The atomic benchmark system: The hydrogen atom H 

First, we test our approach on the hydrogen atom in intense laser fields. For this system 
it is possible to employ a very accurate description of ionization without any absorbing 



boundary conditions using huge basis sets j^J, 

discrete and 2816 continuum states (build from hydrogen eigenfunctions with m — 0) to 
investigate the laser induced dynamics of 



451 ] . E.g., Hansen et al. have used 120 
Id trom hyd 

H(is) y. 



In contrast, we use here a basis of only 40 functions but include absorbing boundary 
conditions. The basis set contains the Is, 2s and 2p z eigenfunctions of hydrogen extended 
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with chains of 37 Gaussians in the direction of the electric field of the laser 54J. The same 



basis has also been used in our previous ca 



however, without absorber. It was found 



culations of H(Ts) in intense laser fields 



[63], 



631 ] that the initial excitation and ionization are 



described very well but fail once the ionization probability becomes noticable (see left part 
of figure HJ) . 

In these and the present calculations, as well as in the reference calculations of Hansen 

n 

et al. |45J the hydrogen atom is exposed to a laser field of the form 

E(t) = E f{t) sin(W + <p) (39) 

with the shape function f(t) 

f sin 2 (£t) for < t < 2T 
f(t)=\ K2Tj , (40) 

otherwise 

T the duration of the laser pulse, a> the frequency, <p the phase and Eq the amplitude. In our 
calculations without absorber jf^] and in the calculations of Hansen et al. [3] the ionization 
probability is defined as the part of the electronic density that is in states with positive 
energy. In our present calculations with absorber the ionization probability is defined as 

P ion = 1 - iV(t final ) (41) 

where N(t^ na {) is the norm of the wave function at the end of the calculation with tg na i = 
2 T + 500 a.u. The additional time interval of 500 a.u. ensures that the norm has definitely 
reached its plateau after the laser pulse. 

In figure H the ionization probabilities of H(ls) as a function of the laser pulse duration 
are shown. The results on the left/right side of figure Q have been obtained in calculations 
without /with absorber potential. Two intensities and two wavelengths are considered. The 
high precision data of Hansen et al. j3] are included. As mentioned already, good agreement 
between our calculations without absorber and those of Hansen et al. is found only if the 
laser pulse is relatively short or weak, i.e. if the ionization probability is small (left graphs 
of figure Q) . The agreement for longer or more intense laser pulses is dramatically improved 
if the absorber is included (right graphs of figure EJ). 

Thus, it is possible to use a small basis set together with absorbing boundary conditions 
instead of an accurate treatment of the continuum. This is of special importance for molec- 
ular many-electron calculations. For these systems it is essential to use small basis sets and, 
thus, to introduce absorbing boundary conditions. 
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without absorber 



with absorber 




T [a.u.] 



>? 10' 




T [a.u.] 



FIG. 1: (Color online) The ionization probability of H(ls) as function of the laser pulse duration 
calculated without (left) and with (right) absorber potential. Two frequencies {oj = 0.55 a.u. 
(top) and u) = 0.18 a.u. (bottom)) and two laser intensities (8.78 x lO 13 -^- (lower curves) and 



3.8 x 10 



15 W 



(upper curves)) have been used. Our results without absorber (Uhlmann et al. 



and with absorber (this work) are compared to the high precision data of Hansen et al. 
numbers in brackets in the legends denote the size of the basis sets. 



63]) 



45J. The 



B. The molecular benchmark system: The aligned hydrogen molecular ion H,j 



Next, the absorber is tested on the molecular benchmark system, pre-aligned Hj, i.e. 
the molecular axis is oriented parallel to the electric field of the laser and, thus, the nuclear 



motion is restricted to this axis. For t 
been performed by Chelkowski et al. 



lis system, full quantum mechanical calculations have 



at ions. 



481 1 which can serve as reference calcu 
In contrast to our previous studies with a minimal 3| and extended |s4 1 LCAO (linear 
combination of atomic orbitals) basis we use here an elaborate basis set consisting of un- 
contracted Gaussians only (details are described in appendix IB 1)1 . It delivers an excellent 
description of the ground state surface with an equilibrium distance of R eq = 1.9975 a.u. 
and a minimum at E min = —0.60246 a.u. (see fig.^J). The ground state energy £"0 = E min + ^ 
and the vibrational levels E n (n = 1,2,...) are also shown in fig. El They are calculated 
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according to the Bohr-Sommerfeld formula 



p(R, E)dR = (n + n ) x h 



(42) 



with uq = 1/2 for the harmonic oscillator, n the vibrational quantum number, h the Plancks 
constant and p(R, E) the classical momentum as function of the energy E and the distance 
R. This yields a binding energy of Eq = —0.59724 a.u., which is of quantum-chemical 
accuracy (cf. e.g. jo^t). 

-0.45^ 



-0.5 



m 



-0.55 



-0.6 







— groundstate surface 
- vibrational states 

- v=6 




FIG. 2: (Color online) Groundstate curve of Hjj" (black) and vibrational levels (thin red lines, thick 
blue line for v = 6) 

The is now exposed to the quasi-cw laser with the parameters according to that of 
Chelkowski et al. 4sj|- So, the laser has a short turn-on of 1 fs. The shape is kept constant 
afterwards. The frequency is to = 0.21 a.u = 5.71 eV and the intensity is 3.5 x 10 13 W/cm 2 . 
To obtain probabilities, 1000 trajectories were calculated and the results were averaged. 
The initial conditions of the trajectories were chosen according to the classical distance 
distribution in the 6th vibrationally excited state. Probabilities are defined as an average 
over the respective quantity for a single trajectory, i.e. 

1 - 

-P<quantity> = ~ ^ ] -P<quantity> ■ (43) 



i=l 
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The ionization probability for one trajectory % is defined as the missing part of the norm, 
i.e. 

PUt) = 1 - Ni(t) . (44) 
The fragmentation probability is defined as 

, for R(t) < R D 
PL g (t)={ (45) 
1 otherwise 



with i?D = 9.5 a.u. taken from |48j]. In accordance with Chelkowski et al. a dissociation 
probability, i.e. fragmentation without ionization, is defined as 

P^ B (t) = (1 - PL®) P L(t) • (46) 

The resulting ionization probability is shown in comparison to the full quantum mechan- 
ical results of Chelkowski et al. j3j in the upper part of figure El The present calculations 
result in an ionization probability that is slightly higher than the full-quantum mechanical 
results. This is due to the different definitions of the ionization probabilities. Whereas in 
grid calculations electronic density can be counted as ionized only at large distances from 
the center of the system (e.g. outside a cylinder with R < 32 a.u. in 0|), our absorber (JHOJ) 
acts also in the vicinity of the nuclei. Therefore, the onset of ionization is somewhat earlier 

n 

in the present calculation as compared to the results of Chelkowski et al. |48[. After 12 fs 
both ionization probabilities have a nearly linear slope. Regarding the uncertainties result- 
ing from the different definitions of absorbing boundary conditions, all in all, very good 
agreement between the present and the full quantum mechanical calculation is found. 

In the lower part of figure El the associated dissociation probabilities are shown. As seen, 
the onset of fragmentation is slightly delayed in the NA-QMD calculations as compared to 
the full quantum mechanical result. This is clearly due to the classical description of the 
nuclei. Afterwards, the same behavior is observed, i.e. first a steep rise and after 18 fs the 
dissociation probability seems to run into a plateau. In spite of this qualitative agreement, 
the present calculation overestimates the dissociation probability by 0.08. 



C. Non-aligned and H2 molecules 

In this section, ionization probabilities and rates for H^" and H 2 are considered as function 
of the angle between the molecular axis and the laser polarization axis as well as as 
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function of the internuclear distance R between the nuclei. Due to the lack of unaligned 
studies, comparison with previous work can be done only for aligned molecules, i.e. = 







49 



50 



70 



71 



72 



73]. 



All calculations have been performed with fixed nuclei. The H 2 molecule is considered in 
the TDHF approximation (see [3] for details). To account the orientation dependence 
the local basis functions centered at the nuclei are extended with functions located at 
hexagonal grid points (see appendix IB 21 for details). With this basis the ground state 

H + 

surfaces exhibit equilibrium internuclear distances and energies of i?eq 2 = 1.99744 a.u., 
= -0.602455 a.u. and = 1.39384 a.u., E^ n = -1.13608 a.u. 
a. Ionization probabilities In the following, the ionization probabilities of and H 2 as 
function of will be discussed. They are calculated at the equlibrium internuclear distance 
with a 50 fs, 266 nm and 5 x 10 14 W/cm 2 laser pulse. For H^, the ionization probability 
Pi on is again defined as the missing part of the norm, eq. (|41j). For H 2 , we abbreviate the 
missing part of the norm of the single-particle functions with 

P s ,i/2 = 1 - N 1/2 (47) 

where N\/ 2 are the norms of the single-particle wave functions of the two electrons. Because 
both electrons differ only by their spin degrees of freedom one has of course N± = N 2 
and P Sj i = P S)2 = P s - The single and double ionization probability are obtained via the 
single-particle approximation [zj, 13] 

Pingic = (1 - NJN2 + iVi(l - N 2 ) , (48) 
Pdoubic = (1-Aq)(l-7V 2 ). (49) 

Note, that with the definition ()48)1 the maximum of P S i n gic is 0.5 if N\ = N 2 . 

P ion , P s , Pgingie and Pioubie are plotted as function of the angle in fig. HJ All calculated 
probabilities have been fitted according to 

n(0) = (P|, cos 2 + Pi sin 2 0) . (50) 

This parameterization has been used by Litvinyuk et al. j3] to fit their experimental data 
on N 2 . The parameterization (|5T)|) fails to fit the experimental data on 2 Q]. As can 
be seen from fig. the parameterization (J50|) works very well for all probabilities in and 
H 2 where experimental data do not exist to date. 
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From fig. 0] it becomes apparent that the probabilities Pi on and P s are largest at par- 
allel orientation = degree and decrease with increasing exhibiting the minimum at 
perpendicular orientation = 90 degree, for both molecules, as expected. However, the 
orientations dependence of the response is much more pronounced for H;j~ as compared to 
H 2 (note the different absolute values of the probabilities in the upper and middle part of 
fig. EJ). In addition, the single ionization probability in H 2 is practically orientation inde- 
pendent, at least for the chosen laser parameters where the double ionization probability is 
not small (lower part of fig. HJ). 

Striking differences between and H 2 have been found also in the alignment behavior 
of both molecules In a relatively large range of laser intensities, fragments originating 



from are much more aligned as compared to those from H 2 76J. This is in accord with 
the present findings of the much larger anisotropic response of in comparison to H 2 . 

b. Ionization rates Finally, the ionization rates T(s~ 1 ) as function of the internuclear 
distance for parallel and perpendicular orientations will be considered for both molecules. 
They are calculated from the logarithmic decrease of the norm N(t), i.e. 

\nN(t) = -Tt. (51) 

In the case of H 2 one has r = Ti + T 2 in the spirit of the independent particle model (i.e. 
N(t) = Ni (t) iV 2 (t) ) . The cw-laser used has a short turn-on of three optical cycles and a 
constant shape afterwards. The wavelength is 266 nm and different intensities have been 
used. For the intensity of 8 x 10 13 W/cm 2 the resulting ionization rates as function of the 
internuclear distance R are shown in fig. El for parallel and perpendicular orientation and 
both molecules. 

For and parallel orientation (left upper part of fig. |5J), the well-known features are 
recovered with our new formalism, i.e., enhanced ionization is observed for internuclear 
distances between 6 and 10 a.u. This is the well-known charge-resonance enhanced ionization 
(CREI) 0Q 

which is accompanied by an electron localization 

0,HQ- Th e CREI 



features are observed here although the system is in the multi photon regime, i.e. the Keldysh 
parameter (see e.g. 4^j]) is 7 ^> 1. This emphasizes the generality of the CREI mechanism. 
At R ~ 3.5 a.u the one-photon resonance between the electronic groundstate and the first 
excited state manifests itself as a peak in the ionization rates in parallel orientation. An 
additional peak is found at R ~ 4.5 a.u. This peak is suppressed for higher intensities, for 
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which the results are not included in figure 03 

For H2 and parallel orientation (left lower part of fig. EJ), the ionization rate exhibits two 
peaks which look similar to the CREI peaks of in parallel orientation. These peaks are 
located at smaller internuclear distances than for Hj. In addition, one of the peaks is split 
which might originate from resonances as the system is in the multiphoton regime. 

The occurrence of enhanced ionization in H2 has been also subject to a number of pre- 
vious studies. First, the distance dependent ionization of H 2 was investigated using a ID 
model jsj. Depending on laser intensity and frequency either one or two peaks were found. 
E.g., at A = 532 nm and I = 1 x 10 14 W/cm 2 enhanced ionization (EI) was found at 
R ~ 4 a.u. and R ~ 6 a.u. which is very similar to the findings in our fully 3D calculations. 
The behavior of parallel aligned H? in static electric fields was studied with high precision 
quantum chemistry methods |tiL 0] , too. An avoided crossing of the H2 groundstate and 
an excited state, that corresponds to the ionic fragments H + and H~ in the static electric 
field, was found. In a fully 3D study of parallel aligned H2 it was found as well 

that the enhanced ionization is linked to the formation of ionic components H + and H~ as 
a typical signature of CREI in the H2 molecule. 

We note also, that the direct experimental observation of CREI in has been published 



y 79j . For H2 this is not the case. Moreover, it has been shown recently in a ID 



recent 

study [73[ that ionization of H 2 usually takes place near the equilibrium internuclear distance. 
Thus, the role of the CREI mechanism in aligned H2 remains still an open problem for future 
investigations. 

The most surprising result of the present studies concerns the ionization probabilities 
in perpendicular orientation (right part of fig. EJ). Clearly seen, the calculations predict en- 
hanced ionization for both molecules. In , the ionization rate exhibits a distinct maximum 
at R ~ 7 a.u., whereas in H 2 , this quantity clearly increases around R ~ 2.5 a.u. (Note also, 
that tiny peaks are found in H2 near the equilibrium distance R eq for both orientations, and 
cf. discussion above). One can definitely exclude the CREI mechanism to be responsible for 
the enhanced ionization in perpendicular orientation, because electron localization cannot 
take place in this geometry. In , this feature probably originates from a resonance that 
occurs around R ~ 7 a.u. However, this as well as the nearly structureless enhancement of 
the ionization rate in H 2 remain as subjects for systematic future studies. 
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IV. CONCLUSIONS 



It was the main aim of this work to present, for the first time, a method to introduce 
absorbing boundary conditions in calculations of many electron systems with basis expan- 
sion. The basic idea is rather general (section 111 AJ) and based on an imaginary potential 
constructed as a projection operator with time-dependent adiabatic states. In this paper it 
has been used to extend the NA-QMD formalism in order to describe realistically ionization 
processes in many electron systems. 

The method was tested on the benchmark systems, the hydrogen atom and the aligned 
H^" molecule in intense laser fields. Very good agreement between the calculated ionization 
probabilities and that of (numerically very extensive) reference calculations has been found. 

The extended NA-QMD formalism allowed us, also for the first time, to study the ioniza- 
tion process in unaligned and H 2 molecules. A completely different orientation depen- 
dence of the ionization probabilities in both molecules was found, with a distinct anisotropic 
response in and a nearly isotropic behavior in H2. In addition, enhanced ionization for 
perpendicular orientation in both molecules is predicted, which to the best of our knowledge, 
has not been reported before. 



The present method can be applied to larger systems, like N2 |32j ]. O2 |33[ or organic 
molecules where the experimental findings still require a consistent interpretation. Finally, 
we note also that the present method can also be used to study ionization processes in 
atomic collisions with many electrons [sol ]. 
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APPENDIX A: BASIS TRANSFORMATIONS 

In order to solve the time-dependent eigenvalue problem (|2T)|) 

H^it)\xV){t) = el(t)\ X :)(t) (Al) 
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the effective single-particle "field-following" adiabatic states are expanded in the basis 

p=i 

Multiplying \<f) a ) with \Xa)(Xa\ the expansion of \<p a ) in the basis {|Xa)l results 

N h 

\<Pa)=Y, S ^a\x:)- (A3) 

a/3 

Note, that \Xa)(Xa\ can b e use d like an identity because the basis sets and {|x a )} 

span exactly the same part of the Hilbert space. Furthermore, the property 

N b 

0=1 



is obtained by using (|A2|) and (|A3|) . Therefore, the transformation £7°" is unitary only if 
both basis sets are orthogonal, i.e. if also (S 1-1 )^ = 5p T 

Inserting (jA2|) into ()A1|) and multiplying with (<f> a \ results in the generalized eigenvalue 
problem 

(MH:^) - e^0 Q |^)J U: p = . (A5) 

(3=1 

The effective single-particle energies e u a and the transformation U a are obtained by solv- 
ing (lA"5l) . 

The effective single-particle wave function l^/- 5 " 7 ) (|15|) can either be expanded in the basis 
{\<p a )} or in the basis {\Xa)} 

|* iff )(t)=E^(*)l^>(*))=E^(*)lx.>(*)) with J = l--.^,ff=U- (A6) 

«=1 a=l 

Using equations (|A2|) and ()A6|) the transformations for the coefficients a£ T and a?° 

< = (x:\*n = E *c< 5 <* > 

a/3 

iVb 

< = E (A8) 

0=1 
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and a general matrix 0°g and O a ah 

O a ab = (Xa\d\ Xb ) = U «aO^ , (A9) 

a/3 

= (<p a \o\<p p ) = Y J s ai u: 1 o: h u^ b + s sp (mo) 

ab-yS 

are obtained. The last transformation is used to calculate the matrix elements V° hs a a which 
are used in practical calculations. 



APPENDIX B: BASIS SETS 

In this appendix all details of the basis sets used in the H^" and H 2 calculations are 
given to make the calculations comprehensible. The basis set used in the H calculations 



(section IlII AJl consists of the hydrogen Is, 2s and_2p z basis functions extended with chains 
of s-type Gaussians and is described in detail in [63] . 

1. Aligned 

The basis set that is used in the calculations to aligned H^" ( sect ion [ill Bj) is a combination 
of a local basis centered at each of the two nuclei and a chain of additional s-type Gaussians 
located along the laser polarization axis. Gaussians 

0Wr') = NY lm (6', ( f>')e ^ (Bl) 



are used also as basis functions at the nuclei. In f)Bl|) Yi m are the spherical harmonics, N is 
a norm constant, <7j, I and m are parameters of the basis functions and f\ = r — Ra where 



Ra is the center of the basis function (see e.g. 



. ail)- 



At each of the two nuclei a basis set 



that consists of such Gaussians is located. The cr, are determined with 

<n = °x f~ l (l<i<N) (B2) 

where a\, f, N are parameters given in table HI With this basis set the atomic as well as 
the molecular groundstate are described extremely well. Furthermore, a good description of 
the excited states of H^" is achieved with this basis set. Please note, that the a% for I = 1 
and I = 2 have been chosen in such a way that a l ^ x = 1.7 1 / 3 <r l =^ x = 1.7 2 / 3 cr l ^ x . 
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f 

1.7 
1.7 
1.7 



<7i [a.u. 
0.05 
0.8473 
1.7191 



fmax [a.U. 

3.487 
4.162 
4.968 



N 
9 
4 
3 



TABLE I: Gaussian basis centered at each of the protons of H^. The parameters given are those 
of equation ()B2|) . 

The parameters of the additional chain of s-type Gaussians (see e.g. that is 

laid out symmetrically to the origin along the z-axis are given in table |H] These additional 
functions have nearly no influence on the already excellent description of the groundstate and 
the lowest excited states. They do, however, improve the description of highly excited and 
ionized electronic states and a dense level structure around E — results. The parameters 
given in table ITT1 were determined using the formalism described in j&J. 



a [a.u.] 


d [a.u.] 


n 


5.54 


3.7 
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TABLE II: Parameters of the chain (width a, spacing d between neighboring functions and number 
n of functions) of s-type Gaussians laid out symmetrically to the origin along the z axis. 



2. Unaligned and H2 



The same basis as before (see table H} is located at each of the hydrogen atoms in the 
calculations to unaligned and H2 ( section fill C|) . However, a hexagonal grid of additional 
basis functions in the y-z plane 

1 ^ 

(B3) 



■ I J 
Zi 



\ % / 





(*-f)f 

^ (j -f)V3d+|(i-f) mod2| J 
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with 



< i < Ni and 
< j < N 2 - 1 if Ni + % even or 
< j < N 2 if Ni + i odd 

is used to extend the basis sets at the nuclei instead of the chain of Gaussians described 
in appendix IB 11 It is thus possible to calculate the response of unaligned molecules to 
the laser field. The parameters of the two hexagonal grids and the one chain are given in 
table UTTl and were also determined using the formalism described in j^]. The last set of 
Gaussians positioned at different places in space has such a large width and therefore also a 
large spacing that it is not necessary to build a hexagonal grid for these Gaussians. Instead 
these basis functions are again laid out chain-like along the z axis |o^ . 



a [a.u.] 


d [a.u.] 


Ni 


N 2 


5.74 


5.2 


9 


7 


7.81 


10.38 


5 


3 


16.62 


18.68 




3 



TABLE III: Parameters of the hexagonal grids (first two lines) and chain (bottom line) of s-type 
Gaussians laid out in the y-z plane. 

It is reasonable to use a basis set for the density in the H2 calculations. The "exact" 
density, i.e. the sum over the absolute square of the single-particle functions, is expanded in 
this density basis to accelerate the calculation of Coulomb and xc matrix elements. Please 
note, that the norm of the density basis functions is different from the usual Lo(R 3 ) norm 

n 

(see [81]). The parameters of the density basis used in the H 2 calculations are given in 
table a x = 0.05/^ has been chosen because the multiplication of two Gaussians with 
the width 0.05 at the same place results in a Gaussian with a width 0.05/-S/2- This density 
basis had not been transformed, i.e. the uncontracted Gaussians were used. Density basis 
functions are also located at the additional centers specified in table IHII The (Jdens of the 
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1 


f 


a\ [a.u.] 


o"max [a.u.] 


N 





1.7 


0.05 
V2 


= 0.035355339 


2.46 


9 



TABLE IV: Gaussian basis used as the H density basis. The parameters given are those of equa- 
tion (lB2l) . 

density basis functions are set to a/y/2. 
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FIG. 3: Ionization (top) and dissociation (bottom) probabilities of as function of time obtained 
with the present NA-QMD method and by Chelkowski et al. . (see text for details) 
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FIG. 4: Ionization probabilities Pi on for ^nd H2 (top and bottom). Also, the missing part 
of the norm of the single-particle wave functions P S) i = P Sj 2 = P s in H2 is shown (middle). The 




FIG. 5: Ionization rates as function of the internuclear distance for (top) and H2 (bottom) at a 
laser wavelength of 266 nm and an intensity of 8 x 10 13 W/cm 2 . The vertical broken line indicates 
the equilibrium internuclear distance. 
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